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Abstract 



In this paper, we propose a method to simulate the microflows with Shakhov model 
using the NPura method developed in [H [6j [5] . The equation under consideration is the 
Boltzmann equation with force terms and the Shakhov model is adopted to achieve the 
correct Prandtl number. As the focus of this paper, we derive a uniform framework 
for different order moment systems on the wall boundary conditions, which is a major 
' difficulty in the moment methods. Numerical examples for both steady and unsteady 

problems are presented to show the convergence in the number of moments. 
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1 Introduction 

(N 

In the kinetic theory, the degree of rarefaction of a gas is often characterized by the 
dimensionless Knudsen number Kn = X/L, where A is the mean free path and L is 
' the relevant characteristic length. The classic Navier-Stokes-Fourier (NSF) equations are 

accurate only when Kn < 0.01. However, the ongoing miniaturization of technical devices 
C^) , requires modelling of gas in microscopic channels, for which the characteristic length L 

is so small that even under normal density and temperature, the Knudsen number is 
beyond the available region of NSF equations. Meanwhile, in the transitional regime 
(0.1 < Kn < 10), the traditional no-slip wall boundary condition is no longer valid. In 
order to match the physical experimentation, the interaction between wall and gas should 
' be carefully conducted. We refer the readers to [13] for more details. 

For microflows, it is known that the Boltzmann equation with Maxwell boundary 
conditions |15| is able to accurately describe the flow state. However, on the computational 
perspective, the cost for solving Boltzmann equation directly is unacceptable in the general 
case. Grad [7] did a pioneer work which extends Euler equations to a thirteen-moment 
system, which opened a new way for modelling rarefied gas flow called as moment method. 
However, it was discovered by Grad himself in [8] that this system fails to give smooth 
shock profiles when the Mach number is larger than 1.65. To remedy this drawback, some 
authors tend to construct a parabolic system similar to the NSF equations. In this field, 
some methods such as Jin-Slemrod [12], COET [IT] and R13 [£l [23] were subsequently 
raised. Concurrently, increasing attention is attracted to systems with more than 13 
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moments (e.g. [23 EI])- As a combination of these two directions, R20 and R26 equations 
were respectively studied in |16] and [11]. In [JJ, a general method for numerically solving 
the regularized moment equations of arbitrary order was proposed, and it was improved 
in [5l [6] and abbreviated as NRxx method in [6] for convenience. On the other hand, 
the boundary condition for the moment methods is a major obstacle for applications of 
moment methods in the field of microflows. In Grad's paper [7J, the basic idea for the 
modelling of Maxwell boundary conditions in the framework of moment method is raised. 
The idea was also used in [28U11] for R13 and R26 equations. However, for general moment 
equations, the numerical method to process the boundary conditions for NRxx method is 
unavailable yet. 

The major concern of this paper is to supply suitable boundary conditions for NRxx 
method. Before that, the NRxx method is first improved such that it is able to predict 
stress and heat flux correctly in the dense case. This is achieved by replacing the BGK 
collision model [2] used in [H0E] with the Shakhov model [19]. Recall that for the BGK 
model, collision term can be analytically solved when using the NRxx method. Similarly, 
analytical solution for each moment can also be obtained when using the Shakhov model. 
At the same time, the force term is also applied to the NRxx method, and one can find 
that this term only affects the momentum equation that it is turned to be trivial when 
splitting method is employed. 

As to the wall boundary conditions, we follow the idea of Grad [7J and try to ap- 
proximate Maxwell boundary condition using moment method. The Maxwell boundary 
condition is a linear combination of specular reflection and diffusive reflection. According 
to Grad's theory, only the moments of odd order in the normal microscopic velocity are 
controlled by boundary conditions. These moments for the specularly reflective part van- 
ish. For the diffusive reflection, the incidence part and the emergence part are considered 
separately. For the incidence part, one need to calculate the moments of a distribution cut 
off by a half space. Since the distribution is expressed by a finite expansion of Hermite 
series, the cut-off turns out to be quite intricate. We eventually derive a simple recursive 
formula to obtain these moments with careful investigation into the detailed expressions. 
The obtained formula brings only slight increment of the computational cost. For the 
emergence part, which is a half Maxwellian, the moments are obtained by direct integra- 
tion, and the result is also given in a recursive form. The overall boundary condition is 
the summation of both the specular part and diffusive part, which is rearranged into a 
simple formulation. It is numerically implemented by first constructing a set of moments 
satisfying the boundary conditions, and then approximating the flow state in the ghost cell 
with a first order extrapolation of each moment. Thus, boundary conditions for the NRxx 
method of all orders are collected into a uniform framework, which avoids separate and 
involved implementation for different systems with sophisticated expressions [26} l28j [TT] . 

A number of numerical examples are presented to show the validity of the boundary 
conditions. Both steady and unsteady problems are studied. Numerical simulations up to 
455-moment system are carried out. The classic symmetric planar Couette flow and force- 
driven Poiseuille flow are investigated as examples for steady problems. All the numerical 
results exhibit the convergence of the NRxx method as the number of moments increases. 

The layout of this paper is as follows: in Section [21 we give a brief introduction to the 
Boltzmann equation and the NRxx method. In Section [3l the Shakhov collision model 
and the force-induced acceleration terms are coupled with the NRxx method. In Section 
HI the derivation of boundary conditions are carried out. Numerical examples are shown 
in Section O and some discussions on the validity and accuracy of the NRxx method are 
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given in Section [6l Finally, we make some conclusion in Section \7\ 



2 The Boltzmann equation and the NRxa? method 



The Boltzmann equation is the basic equation in the kinetic theory, where a distri- 
bution function f(t,x,£) is introduced to provide a statistical description for the motion 
of molecules. Here t £ K + is the time, and x, £ 6 M 3 are the position and velocity of 
particles. The Boltzmann equation reads 



^ + S-V x f + F-V $ f = Q(f,f), 



(2.1) 



where F is the acceleration of particles caused by external forces. The detailed expression 
of the collision term Q(f, /) is not presented here due to its complexity, but we stress 
that Q(f,f) contains a five-dimensional integration which causes great difficulty in the 
numerical simulation. Instead, simplified collision models such as the BGK model [2] and 
the Shakhov model |19] are adopted in this paper. These models read: 



1. BGK model: 



2. Shakhov model: 



dt 



1 



+ £-V !B / + F-V£/ = -(/m-/); 



(2.2) 



% + t-V xf + F.Vsf= 1 - 



1 + 



(l-Pr)(£-u)- q (\i 



5 P 6 2 



u 



Jm - f 
(2.3) 



Here p, u, 8 and q denote the density, mean velocity, temperature and heat flux respec- 
tively, and these macroscopic variables are related with the distribution function / by 



P 



u\ 2 fdt 



1 

u = - 

p 

1 

q= 2 



3 

z-u\\z-u)fda. 



(2.4) 



Besides, r is the relaxation time and fu is the local Maxwellian which can be analytically 
formulated by 



Jm 



P 



exp 



u 



29 



(2.5) 



(27r0) 3 / 2 

In (|2.3p . Pr stands for the Prandtl number which is a constant. One can easily observe 
that if Pr = 1, then the Shakhov model reduces to the BGK model, which agrees with the 
common knowledge that the BGK model predicts an incorrect Prandtl number 1. 

The NRxx method is a numerical tool for solving large moment equations. It originated 
in [1] and was simplified in [5]- The basic idea is to expand the distribution function / 
into the Hermite series: 



f(t,X,£) = fa(t,x)n e ,a I 

aeN 3 V 



£ — u(t, x) 



(2.6) 
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where Hd t a is the basis function defined as 

3 



1 / 2 \ 

d=i v 27r V 1 J 



(2.7) 



and He n is the Hermite polynomials 



He n {x) = {-IT exp (^) exp (-£) . (2.8) 

For convenience, we let He n {x) = if n < 0. Thus He^ a {v) is zero when any of the 
components of a is negative. 

With the expansion (|2,6p . the coefficients / Q can be considered as a set of infinite 
moments, and we have the following relations: 



fo = P, fe z =0, he d = 0, 

d=1 3 (2.9) 

&ij = fet+ej, &ii = 2/2^, ft = 2/3 ei + f2e d +e t , 



d=l 



where i, j = 1, 2, 3 and i 7^ j, and ajj is the stress tensor or pressure deviators, which can 
be deduced from the distribution function / by 

(Tij = Pij - -Sij^pda, with Pij = [ (£i-Ui)(£j -Uj)fd£, i,j = 1,2,3. (2.10) 
3 d=i J * 3 

In order to implement (|2.6p numerically, a positive integer M ^ 3 is chosen and only 
the coefficients {f a (t,x)}\ a \^ M are stored. Due to the absence of higher order moments, 
the resulting moment system is not closed. According to [5], the (M + l)-st order moments 
are approximated by 



- T i p 2^ dx ^ + D \2^ dx 2^ 2^ 

l' j=l J \j=l 3 J d=l j=l 

D 



9f a - 



dxj 



E( dud 1 d9 ... , . , 

V ~dx~ * a ~ ed ~ £j + 2dx~^ * a ~ 2ed ~ ej + + l >J a - 2e i+ e i 
d=i ^ J J 



(2.11) 



Here / Q is taken as zero when any of a's components is negative. The numerical scheme for 
the force-free BGK model has been constructed in [6j based on the finite volume scheme 
with linear reconstruction and the fractional step method. Suppose the problem is in ID 
and the grid is uniform with cell size Ax. We denote the cell centers as Xj, and then a 
full time step of the scheme can be sketched as follows: 

1. Determine the time step size At. 

2. Reconstruct the first M-th order moments for the distribution functions on cell 
boundaries Xj±i/2 with a conservative linear reconstruction. 

3. Get the (M + l)-st order moments for the distribution functions on cell boundaries 
with a direct discretization of (|2.1ip . 
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4. Apply the HLL scheme to solve the purely advective equation dtf + £• V^/ = over 
a time step of length At. 

5. Analytically solve the pure collision equation of the BGK model dtf = (/m — /)/ T 
over a time step of length At. 

We refer the readers to [2J [5j [6] for details of the algorithm. Here we only note that the 
Step [5] is nontrivial since two distributions cannot be added up directly, and in Step El the 
reason why the collision-only equation can be directly solved is that /m can be expressed 

in the Hermite series {%e,a\ trivially as fu = faHo,a — u)/yf~0\. 

3 The NRcccc method for Shakhov model with force terms 

As is well known, the Prandtl number for monatomic gases is around 2/3, while the 
BGK model gives a Prandtl number 1, which causes incorrect prediction of the stress 
tensor or heat flux q for a dense gas. As a remedy, the Shakhov model was introduced 
in |19] as a generalization of the BGK model. The difference between these two models 
has been investigated in [SHE]- In this section, we extend the NRxx method in [5] to 
the Shakhov model, and the force terms in (j2.3[) is added. 



3.1 The governing equations 



The moment system for the Shakhov model (|2.3|) with moment set {f a (t, a;)}|a|^M wm 
be deduced here. As in [5], the strategy is to expand (|2.3p into Hermite series, and then 
match the coefficients for the same basis functions. In order to simplify the notation, we 
define 



df 

B = F ■ V € /, 



(3.1) 



C= \ { 



[I - Pr)(£ - u) ■ q (\£-u 



5p9 2 



fu - / 



It has been deduced in [5] that the Hermite expansion of A is 

\ d=l d=l / 



aeN 3 

3 



i=i 



dfa- 



d Xj 



9fa- 



d Xj 



du 

+ Y i e fa-e d - ej + Ujf a -e d + (pLj + l)f a - ed+ej ) 

d=l i 



1 80 



+ 2~dx~- Y { e fot-2e d -ej + Ujf a -2e d + (otj + l)f a ^ 2 e d +e 



d=l 



Using the differential relation of the Hermite polynomials, we have 



_d 



ti-u 



,a+e d 



Ve 



(3.2) 
(3.3) 
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Thus the Hermite expansion of the force term B can be easily deduced as 

3 



B ~ ~ Yl ^2 F df a -e d n eia (—7=-) 



(3.4) 



The expansion of the collision term C can also be obtained by direct calculation. The 
result is 



C 



1 



1 -Pr 



3 3 



i=l j=l 



^ faH-e,, 

\a\^2 



sTq 



(3.5) 



Putting (|3.2p(|3.4p and f)3 . 5 j) into the Boltzmann-Shakhov equation A+B = C and extract- 
ing coefficients for all basis functions, with a slight rearrangement, we get the following 
general moment equations for Shakhov model: 



9f a | uu d ^ uu d I i ( ua v-^ w 



<9u„ 



d=l 
3 



St 



dx 



1 / d# 



30 



2 <9i 



^ (#/a-e d - ej + («j + l)/a-e d+ej ) + ( e f*-2e d - ej + (otj + l)/ a _2e d +ej 



+ E 



where 5ij(a) and 5(a) are defined by 



3 



1 / 1 - Pr 



^2 5 ij (a)q i - 5(a) f a 



1, ifa = ej + 2ej, 
0, otherwise, 



5(a) 



1, if|a|>2, 
0, otherwise. 



(3.6) 



(3.7) 



Now we will explore something more from (|3.6|) . Noting that / e . =0, Vj = 1,2,3, the 
following relation can be obtained if we put a = into (|3.6p : 



dfo 



3=1 



ax, 



(3.8) 



This is the mass conservation law. If we set a = e d , d = 1, 2, 3, the equations are 



/o 



dx d dx d 



e d +ej 



3=1 



d Xj 



0. 



(3.9) 



This equation can be simplified as 



3=1 J 



0. 



(3.10) 



Now we consider the case of \a\ ^ 2. Substituting (|3.10p into (|3.6p . the temporal differen- 
tiation of u can be eliminated. In order to eliminate the temporal differentiation of 6, we 
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multiply (|2,3p by |£ — u\ 2 on both sides and then integrate on M 3 with respect to The 
result is 



3 o„\ „3/ 3 



^I + E»^j + |E^ + E^j-. (,n, 

Note that the force term does not appear in this equation, since 

/ |£- M | 2 |fd£ = -2 / (0-^)/d£ = 0. (3.12) 

7R3 O^j J R 3 

Thus, the final form of equations for \a\ ^ 2 reads 
^fa _ dPjd f J.VI^i-LV Ir^M I f 



{Ofa~e d - ej + (aj + l)/a-e d +e 3 ) + {9f a -. 2 e d -e 3 + {otj + l)/ Q _2e d + ej - ) 

A ( df a - ej df a df a+ej \ i / i_p r 3 

(3.13) 

In order to get a closed system, we collect flZHD , IpTIO]) . (jHTTT]) and (f3~T3j) with 2 < |a| < M 
together. Then the governing system for the NRxx method with Shakhov model and force 
terms is formed. 

Remark 1. In the Shakhov model, the equation (|2.11|) . the prediction of f a with \a\ = M+l 
derived for the BGK model, is still available. In [5], (|2.1ip is deduced in the following two 
steps: 

1. Determine the orders of magnitude for all f a using Maxwellian iteration. 

2. For | a | = M + 1, remove all the high order terms in the equations containing only 
—fa/r in their right hand sides. 

The Maxwellian iteration can also be applied to (|3.13|) . and after the first iteration step, 
we immediately get 

f a = 0(r), \a\=2, or a = ei + 2ej, i,j = 1,2,3, (3.14) 

and other moments with \a\ ^2 remain to be zero. This result is the same as that we 
have derived in the BGK model. We note that for \a\ > 3 and a = (1, 1, 1), (|3.13p is just 
the corresponding equation for the BGK model. Thus, in the view of order of magnitude, 
the subsequent iterations are identical to the BGK case. Moreover, when M 3, which 
we have assumed in the last section, step is also identical for both models. Hence (|2.11|) 
still applies for the Shakhov model. 

3.2 The numerical approach 

The acceleration F only appears in (|3.10p in the governing system, thus a splitting 
method can be applied as follows: 
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1. Transportation: solve the force-free Shakhov equation over a time step of length At. 

2. Acceleration: solve dtu = F over a time step of length At. 

In order to solve the force-free Shakhov equation, another splitting of the convection and 
collision part is needed. For the convection part, the method is identical to that used in 
the BGK model. We refer the readers to [H [6] for details. For the collision part, since a 
new collision model is adopted, the procedure is slightly different. 

Now we consider the pure collision model, where p, u and 6 are not changed while time 
evolves. Therefore, the collision terms only exist in (|3.13p with \a\ ^2. Two cases are 
considered below: 

(1) a = ei + 2ej, i,j = 1, 2, 3. In these cases, the pure collision equations are written 



MS 



dt 

(3.15) 



1 

T 




2f3ei + ^ fei+2e 3 ) ~ f ei +2e 
3=1 



i,j = 1,2,3. 



In the general case, r only depends on p and 9. Thus it is invariant in the collision-only 
system. This turns ()3.15p into a linear ordinary differential system with 9 equations, which 
can be analytically integrated as 

f ei +2 ej (t) = jjr<?i(to)exp ^- Pr( ~^_ - Qgi(*o) - fei+2 ej (to)^ exp (^-^y^j , (3.16) 

where to denotes the initial time. 

(2) Other cases. For other a's, the collision-only equation is the same as the BGK 
model: 

it = -;'«■ (317) 

The solution is 

/«(<) = /«(*o)exp(-^V (3.18) 



When (|3.16p and (|3.18p are used in the numerical scheme, we replace t and to with t n+ i 
and t n respectively. Note that when r is independent of u, the acceleration and collision 
do not coupled with each other, thus the splitting is applied only once rather than twice. 
This makes it more efficient when the Strang splitting is employed. 



4 Boundary conditions 

In the moment methods, the boundary condition is always a complicated issue when 
simulating microflows. As discussed in [71 [20j [261 [281 [11] an d the references therein, delicate 
derivations and careful numerical techniques are needed for a solid wall. In this section, a 
numerical way for dealing with boundary conditions in the NRxx method is introduced, 
which appears to be uniform for all orders of moment systems. 
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4.1 The kinetic boundary condition 

In the kinetic theory, the most extensively used boundary condition is the one proposed 
by Maxwell in [15]. According to the common hyperbolic theory, for (j2.3|) . the boundary 
condition is only needed when £ • n < 0, where n is the outer normal vector of the 
boundary. For a point x on the wall, supposing the velocity and temperature of the wall 
to be u w (t, x) and 9 W (t, x) at time t, Maxwell proposed the following boundary condition: 

,u ,> \ X/$(t,x,() + (l-x)f(t,*,e), C w n<0, 
/ft "' <) = \/tt-,e, C»-n>0, <41) 

where x £ [0, 1] is a parameter for different gases and walls, and 

= £ - 2{C W ■ n)n, C w = £-u w (t,x), (4.2) 

w p w (t,x) ( \tj-u w (t,x)\* \ 
m (t, x, $,) = — -— exp — T^- — r — . (4.3) 



(2irO w (t,x))V 2 *\ 29 w (t,x) 

The functions u w (t,x) and 6 w (t,x) are prescribed and stand for the wall velocity and 
temperature at time t and position x, and p (t,x) ensures the conservation of the mass 
at the wall, that is, 



/ {C w -n)f(t,x,$)d£ 
X ( [ (C w - n)fZ{t, x, d£ + / (C w ■ n)f(t, x, £) dA = 0. 

\Jc w n<0 Jc w n>0 ) 



(4.4) 



For this boundary condition, the normal velocity of gas on the boundary is the same 
as the normal velocity of the wall. However, in the case of shear flow, velocity slip and 
temperature jump will appear on the boundary. 

4.2 The boundary conditions for the NPLex method 

The boundary condition can be derived by taking moments on both sides on (|4.ip . 
Before that, we define 

C e , a = K , , , , W>0, a G N . (4.5) 
This definition leads to 

g a = C 9ia [ g{£)U eya {v)^ V {\v\ 2 /2)& Vl (4.6) 

where v = (£ — u)/y/6 and g(£) is a distribution function expanded into Hermite series as 
= Sqgn 3 9a7~Ld,a( v )- In order to simplify the calculation, we suppose n = (0, 1, 0) T . 
Thus, taking moments for (|4.ip requires half-space integration 

Ce, a f g(meA v ) e M\v\ 2 /2)dv. (4.7) 

Suppose an M-th order system is used in the NRxx method; that is, an (M+l)-st order 
approximation of the distribution can be obtained through (|2.1ip . This approximation is 
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directly used in (j4.7|) so that the integral can be actually worked out. Concretely speaking, 
(|4.7p is approximated as 

V gpCe. a I U e:a {v)U e ,p{v)eK V {\v\ 2 /2)Av. (4.8) 



\/3\^M+l 



Since u 2 = on the boundary, the region of integration can be written as {v 2 ^ 0}. 
Thus, we only need to calculate 

IaA d ) = Ce,a[ U e , a (v)HeA v )eM\v\ 2 /2)dv. (4.9) 

The details can be found in Appendix [Bl and the result is 



and 



I a ,p{6) = S(a 2 ,p 2 )e^ ■ 6 aih 5 a3 p 3 (4.10) 



1/2, m = n = 0, 

K(l,n — 1), m = and n / 0, 

S(m,n) = 4 ) m j n ( 4 - n ) 

K{m, 0), m 7^ and n = 0, 

K(m, n) + S(m — 1, n — 1) • n/m, otherwise, 
where 

K(m,n)= K ' He m ^{0)He n (0). (4.12) 
m! 

The above deduction leads to the following proposition: 

Proposition 1. Suppose g(v) is a function defined on M 3 which can be denoted by a finite 
expansion of Hermite basis functions 

g(v)= 9cflie, a (v). (4.13) 

\a\^M+l 

for some 9 > 0. Let g(v) be a half-space cut-off of g(v) as 
Then g can also be expanded into Hermite series as 
where 1^(6) is defined in flUBD — dUZD • 

Proof. It is already known in [3] that {^6»,a(' u )}aeN 3 * s an orthogonal basis of the weighted 
L 2 space L 2 (R 3 ; exp(H 2 /2) dv). Since 



/ |3( W )| 2 exp(|^| 2 /2)d- y = / | 5 (u)| 2 exp(M 2 /2)ch; 
< / \g(v)\ 2 exp(\v\ 2 /2)dv = V C^| 5a | 2 <+oc, 



(4.16) 



|a|<M+l 

g(u) also lies in L 2 (IR 3 ; exp(|?;| 2 /2) dv). Thus the validity of (|4. 15|) can be naturally 
obtained. □ 
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The following proposition depicts the sparsity of I a ,p- 

Proposition 2. // l a ,p{0) is nonzero, then (1) a\ = (2) a-s = (3%; (3) a?2 — P2 is zero 
or odd. When a = f3, I a ,p{9) is equal to 1/2. 

Proof. If l a> p(0) is nonzero, (|4,10p directly gives a± = j3\ and Q3 = j3%. If 02 — @2 is a 
nonzero even integer, K{o.2^2) is zero since He n (0) is zero when n is odd. In order to 
prove I a ^{9) = in this case, according to ()4.10p . we only need to prove 5(a2,/?2) = 0. 
This can be done by induction: 

(1) If a2 = or /?2 = 0, 0.2 and $2 must be both even but one of them must be positive. 
Equation (|4,lip shows 5(02,^2) = directly. 

(2) Suppose S(a.2 — 1, $2 — 1) = 0. Then, according to the last case in (|4.11|) . 5(02,^2) 
is also zero. 

Finally, when a = (3, (|4.1Up gives I a> p(9) = S (02,^2)- The subsequent proof can also 
be done by induction, since 5(0,0) = 1/2 and K(n,n) = for n > 0. □ 

According to Proposition [21 we find that only ([«2/2] + 1) terms are nonzero in the 
summation (|4.8p . This greatly reduces the computational cost. 

Now let us return to the boundary conditions. According to Grad's theory [7J, 
in order to ensure the continuity of boundary conditions when x ~ * 0i only a subset 
of moments {/ Q | \a\ ^ M + 1 and 02 is odd} should be used to formulate boundary 
conditions. This will be completed in the following three subsections. Later in this section, 
for conciseness, the variables t and x are omitted in our statement if not specified, and all 
spatially dependent functions are considered to be on the boundary. 

4.2.1 Determination of p w 

For simplicity, we factorize the right hand side of (|4.ip into three parts and consider 
each part independently. Define 

m = r m\ 6 < <, m = r m 6 ; <, = + 



0, " v>/ lo, 

(4.17) 

Then (|4.ip can be rewritten as 

/(£) = + + (1 " (4-18) 

Suppose the Hermite expansion of / is 

/«)= E /«Wff,af^y ( 4 - 19 ) 
|a|<W+l V V# / 

Then q(£) can also be expanded into Hermite series according to Proposition [1] and [2] as 

<Z(0 = £ fe «,, a . (4.20) 

Substituting ([P]) and (|4.20[> into (|4~3|) . p w can be worked out as 

= V W E 5'(l>2fc)0 1 /2-fc /2fee2) (4.21) 



k=0 



where the expression of q e2 is derived from ()4.10p . (|4.15p and Proposition [2l 
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4.2.2 The moments of p and r 

Now the moments for q(£) have been calculated in (|4.20p . we still need to get Hermite 
expansions of p(£) and r(£). We suppose that p(£) can be expanded under the basis 

Then, according to (|4.3|) . the coefficients can be formulated by 

f p w / \£-u w \ 2 \ f\v\ 2 \ 
Pa = ^ L<o (2^)3/2 eXP V 2^^ J 6XP V~ ) ^ (423) 

where £ = yfOv + Define 



1 »n /-+ 00 1 / |v^y 



X 



2 



JJ X ) = -R— / exp - ' "... ' He s (y)dy, (4.24) 



hj+i f° 1 / |v^y 



,r 



2 



JJ X ) = -0— / exp - ' " ' He s (y)dy. (4.25) 

Then p a can be expressed by 

Pa = p W J ai (uY ~ ui)J a2 (uY - u 2 )J a . A (uf - u 3 ). (4.26) 
J s (x) and J s (x) can be calculated recursively as 

) w - e)J s - 2 (x) + xJs-^x)] , Ol; (4.27) 

? w -0)J s _ 2 (x)+aJ s _i(aO] -F s (x), O 1; (4.28) 

^-8H s ^ 2 (x), s > 2. (4.29) 

- 1) 

The starting values are 

J_i(s) = 0, J (a?) = l, (4-30) 

J_i(x) = 0, J (x) = ierfc (-/==)> ( 4 -31) 



J s (x) 


-7l< 






J s (x) 






s L 


H s (x) 





6> w / 



i/ (x) = 0, Hl (x) = ] j — exp^-^j, (4.32) 

The detailed derivation of (|4.27|) - (|4.32j) can be found in the Appendix [Cj Noting that 
u 2 = , (|4.26|) can be further simplified as 

Pa = p W J ai (uY - ui)J a2 J a3 (uY - u 3 ), (4.33) 



where 



J S = -(9 W -6)J S _ 2 -H S , a>l, H s = - t \ 0H S . 2 , s > 2, 
s sis — 1) 

r4- 4.34 

J^ 1 = H = 0, J = 1/2, ffi = y— . 
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Here we emphasize that due to equation (|4.2ip . all p a : s are only related with {/2fce 2 }o^fc^rAf/2l 
besides u, u w , 6 and 6 W . 

Now we turn to the moments of r(£). Note that only the moments with odd 0.2 
are needed. However, r(£) is an even function with respect to C 2 , which causes all its 
moments with odd Q2 vanished. This indicates that r(£) can be simply neglected when 
discussing the boundary conditions. 



4.2.3 Construction of boundary conditions 

Now we take moments with odd «2 on both sides of (|4.18p . Making use of Proposition 
[2j we have 

fa = XPa + XQa = XPa + ^Xfa + X E ' ^^"^"^./a+^fc-^H , ( 4 -35) 

k=0 



where K 2 (a) = \(M — a\ — 0:3 )/2] . A simple rearrangement gives 

2X 



2 - X 



K 2 (a) 

Pa+ S ( a 2,2k)e^ 2 - k f a+{2k ^ 2)e2 
k=0 



(4.36) 



Equations (|4.36p with \a\ ^ M+l and odd a 2 , together with u 2 = form the boundary 
conditions of the dynamic moment equations. Recalling 

p a =p a {u,U W ,e,e W ,fo,f2e 2 ,--- , h \M/2~\ e 2 ) , (4-37) 

one can find that the terms which appear on the left hand side of (|4.36p never appear on 
its right hand side. Thus, if an arbitrary distribution function denoted as (|4.19p is given, 
we can define a functional F b which maps (glgD to another distribution / 6 (£): 

A£)= E (4-38) 

where u b = (u\, u 2 , U3), 9 b = 6, and 

rb _ f fa, if a 2 is even, ^ _ 

01 1 the right hand side of (|4.36p . if a 2 is odd. 

Thus f a satisfies the boundary condition. The mapping F b will be used in the numerical 
implementation of boundary conditions. 

At the end of this section, we prove that u b and 6 b are the corresponding velocity 
and temperature of the distribution function f b ($). This is equivalent to the following 
proposition: 

Proposition 3. If a distribution /(£) with expression (|4.19p satisfies 

3 

f ei =fe 2 =fe 3 =^2f2e d =0, (4.40) 
d=l 

then f b = F b (f) with expression (|4.38p also satisfies 

/eW£=/e3 = E/le d =0. (4.41) 
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Proof. Equation (|4.39p gives 

f b ei =f ei , f%=fe 3 , fL d = f2e d , d= 1,2,3. (4.42) 

Thus it only remains to prove /g 2 = 0. According to (|4.30p . (|4.33|) and (j4.34j) . p e2 can 
actually be expressed by 



p e2 = P w HuY - U i)JiJo(uf - u 3 ) = p w [(6 w - e)J-x - H X ] = -p W ]J (4.43) 

Since ^2(^2) = \M/2~\, the above equation together with (|4.2ip and (|4.36p immediately 
gives fl =0. □ 

4.3 Numerical implementation of boundary conditions 

In a finite volume scheme, the boundary conditions is often applied by ghost cell 
techniques. Suppose the distribution function of the cell on the boundary is denoted as 
(|4.19p . The distribution function of the ghost cell can be constructed as follows: 

1. Apply F b on /(£) and suppose the result is (|4.38[) : 

2. Construct the ghost cell distribution as 

/ shost (£)= E ( 2 fa-f«)ne,*( *~ {2U *~ U) ). (4.44) 

Now we consider the time complexity of this operation. Suppose Nm = (M + 2)(M + 
3)(M + 4)/6 is the number of moments involved in the boundary condition. Obviously, 
(|4,44p requires O(Nm) operations. For the calculation of F b (f), we list the cost as follows: 

1. Half-space cut-off of / (|4T20|) : 0{MN M ) operations; 

2. Calculation of p w (|4T2Tj) : 0(1) operations; 

3. Calculation of p a (j4.33p : 0{Nm) operations; 

4. Evaluation of (|4.39j) : O(Nm) operations. 

Thus, the total computational cost is O(MNm), while the time complexity is O(Nm) if 
no boundary condition is considered. However, since this procedure only takes place on 
the boundary, it produces little increment of the computational time in real computation. 

Remark 2 . Proposition [3] indicates the conservation of mass on the boundary when using 
the HLL numerical flux as in [6]. One can find that when = 0, saying a special 
reference coordinate system is used, the minimum and maximum signal speeds in need of 
the HLL flux are opposite numbers. Together with p§ host = p } u| host = — « 2 , the mass 
conservation of the HLL scheme follows naturally. 

5 Numerical examples 

In this section, three numerical examples are presented to validate our algorithm. In 
all these examples, a hard sphere gas is assumed, for which the relaxation time is defined 

as 

16 V 9 p ( ' 
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following [3], where Kn is the Knudsen number. The CFL number is always 0.95. And for 
all the tests, the wall is set to be a fully diffusive one (x = 1) with 6 W = 1. The POSIX 
multithreading technique is utilized in our simulation, and at most 8 CPU cores are used. 

5.1 The beginning of a shock wave's formation 

The first example is a simulation of the interaction of a coming flow with a diffusive 
wall. The computational domain is [—5,0] and the global Knudsen number Kn used in 
(|5.ip is set to be 0.5. The left boundary is a free boundary, and the right is a stationary 
diffusive wall parallel to the xz-plane. The initial condition is given by 

PQ (y) = 1.0, u (y) = (0, 0.5, 0) T , 9 (y) = 1.0, Vy G [-5, 0], (5.2) 

and the gas is in equilibrium everywhere. A left-going shock wave will form after a suf- 
ficiently long time. Here we stop the computation at t = 1.0 in order to check the 
validity of the boundary condition. For a reference solution, we solve the Shakhov equa- 
tion (|2.3p directly using a Conservative Discrete Velocity Method (CDVM) introduced in 
|27j . For the computation of both NRxx method and CDVM, a uniform mesh with 500 
grids are used to discretize the domain. For CDVM, the computational velocity domain 
is [-10, 10] x [-10, 10] x [-10, 10] and discretized by 50 x 100 x 50 grids. 

Figure [1] and [2] are the results for CDVM and NRxx method for M = 3 to 12. Only 
the part y G [—3,0] is shown since all variables for the remaining part are almost constant. 
Since a large Knudsen number is considered, predictions from lower order moment equa- 
tions give very large deviations, so the necessity of high order moment theory is obvious. 
As the number of moments increases, all profiles get closer and closer to the results of 
CDVM. When M reaches 11, the density and temperature plots agree with the CDVM 
results very well, and the errors in 022 arid 52 are much smaller than the low order cases, 
though it is still observable. It is reasonable that higher order moments converge more 
slowly than lower order moments, which is also observed in |llj . 

5.2 Planar Couette flow 

The planar Couette flow is a classic benchmark test in the field of microflows. The 
moment method for this problem has been investigated in a lot of papers such as [25} \TE[ 
l26l [T0l [28l [TT] . Here we consider the symmetric Couette flow. The gas lies between two 
plates parallel to the xz-plane. Two plates move in the opposite direction with constant 
velocities within their own planes. A steady state can be obtained for a fully developed 
flow. 

In this example, the computational domain is [—0.5, 0.5]. The velocities of the left and 
right plates are 

uf = (-0.6296, 0, 0) T , u% = (0.6296, 0, 0) T . (5.3) 
The initial state is a global equilibrium with 

Po(y) = 1, u (y) = 0, 6 (y) = 1, My G [-0.5, 0.5]. (5.4) 

The steady state can be achieved if the computational time is sufficiently long. Also, 
both the NRxx method and CDVM are applied to this problem. Three different Knudsen 
numbers, Kn = 0.1,0.5,1.0, are investigated. For CDVM, the computational velocity 
domain is chosen as [—10, 10] x [—10, 10] x [—10, 10], and 50 x 50 x 50 grids are used. Here 
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Figure 1: Density and temperature plots for the problem in section I5TT1 The left axis is 
for the dashed lines, and the right axis is for the solid lines (to be continued). 



we note that such discretization may not produce numerical results accurate enough as 
the reference solution, but the computation is already extremely slow. 

Numerical results for Kn = 0.1 are shown in Figure [3] and [H In this case, most 
lines agree with each other. The convergence in the number of moments can be observed, 
however, due to the numerical error from both NRxx method and CDVM, small deviations 
between the CDVM results and the possible limit of the NRxx method can be found. 
One can disclose that lower order NRxx results deviate from the CDVM results more 
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(i) M = 11 (j) M = 12 



Figure 1: Density and temperature plots for the problem in section I5TT1 The left axis is 
for the dashed lines, and the right axis is for the solid lines. 

than higher order ones. This correctly reflects the behavior of NRxx method under low 
Knudsen numbers, as is also found by [I]. 

Now a larger Knudsen number Kn = 0.5 is considered, and the results are given in 
Figure [5] and [6l In this case, the results for odd and even orders evidently break into two 
groups, and they approach closer to the CDVM results separately. This can be also find 
in Figure Q] and [2l For p and 8, the even group gives better results, while for a\i and 
CT22, t ne odd group is more accurate. The reason remains to be further explored. The two 
subfigures in Figure [6] clearly exhibit the convergence. In [28} lllj. it was discovered that 
the normal stress 022 is difficult to match by R13 and R26 equations. Here one may find 
that when the number of moments is increasing, the quality of the approximation to this 
quantity is improved continuously. In case of M = 9, the profile agrees with the CDVM 
result quite well, and when M = 10, the relative difference is below 5%. 

The severe case Kn = 1.0 is also studied. Similar results with the case Kn = 0.5 are 
obtained in Figure [7] and while the magnitude of the difference is much larger. For 022, 
now the relative difference for M = 9 is about 10%. But the rate of convergence is still 
encouraging — compared with the result with M = 4, the error is halved. 
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Figure 2: Stress and heat flux plots for the problem in section [570 The left axis is for the 
dashed lines, and the right axis is for the solid lines (to be continued). 

5.3 Force-driven Poiseuille flow 

This is another example which is frequently used to verify the boundary conditions 
of moment methods |28t 111] , Similar with the Couette flow, the gas also lies between 
two parallel plates, but the plates are stationary and an external constant force parallel 
to the plates causes the flow to reach a non-stationary steady state. In our settings, the 
computational domain is again [—0.5,0.5], and the Knudsen number is set to be 0.1. The 
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(i) M = 11 (j) M = 12 

Figure 2: Stress and heat flux plots for the problem in section [5TT1 The left axis is for the 
dashed lines, and the right axis is for the solid lines. 



force introduces an acceleration 

F = (0.2555, 0,0) T . (5.5) 

The initial condition is the same as the Couette flow. These settings are the non- 
dimensional form of the test in [32], where the DSMC result is carried out, and this 
example is also considered in \60\ [TT] . Since it is quite difficult for us to exert the force 
term in CDVM, we have to use the DSMC result in [32] for comparison in spite of the 
difference in the collision model. 

The numerical results are presented in Figure [9] and [TUJ For all the profiles, the 
convergence in the number of moments is legible, while the NRxx results do not converge 
to the results of DSMC. This may due to the difference between the collision terms of 



Shakhov model and DSMC. Taking the temperature plot (Figure 9(c) ) as an example, 
the result of M = 3 matches DSMC result best, since when M = 3, the collision term of 
Shakhov model is almost the same as that of DSMC. While when the number of moments 
increases, the collision term deviates away from DSMC's gradually. Here the accuracy of 
collision models is not the topic of this paper. Even though, two results are very close 
quantitatively, which indicates the correctness of the boundary conditions and the Prandtl 
number of the NRxx method. 
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M = 3 

M = 4 




1 .02 1 ' ' 1 1 1 ' ' 1 1 

-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 

(b) Temperature, 6 

Figure 3: Density and temperature plots for the planar Couette flow with Kn = 0.1 (to 
be continued) 

6 Some discussions on the NRa?a? method 

6.1 Order of accuracy 

For the macroscopic equations, a basic quantity describing its ability is the order of 
accuracy with respect the Knudsen number. The definition of order of accuracy can be 
found in textbook |22j : 

A set of equations is said to be accurate of order \q, when the pressure deviator 
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Figure 4: Shear and normal stress plots for the planar Couette flow with Kn = 0.1 



cry and the heat flux qi are known within the order 0{e °). 

Here e is a small parameter proportional to the relaxation time r. As we have discussed 
in Remark [TJ in the view of order of magnitude, the process of Maxwellian iteration for 
the Shakhov model is identical to the BGK model. Hence, for an arbitrary M ^ 3, the 
leading order term of f a with \a\ = M + 1 is known from the corresponding moment 
equations (see [S] for details). And it has been deduced in [S] that f a ~ 0(7" n^l /3l ^ f or a ll 
|a| ^ 4. Thus, from the analytical form of the moment equations (|3. 13|) . we immediately 
have that f a with |a| = M is known up to (|~(M + l)/3] + l)-th order. Subsequently, f a 
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Figure 5: Density and temperature plots for the planar Couette flow with Kn = 0.5 



with \a\ = M — 1 is know up to (|~(M + l)/3] + 2)-th order, and this can be done until 
| a | =2. The general result is 

Proposition 4. For the moment equations described in Section \3.1{ f a has ( [~4(M + 
l)/3] — order accuracy if 2 ^ \a\ ^ M. 



Now, using (|2.9p and the definition of order of accuracy, we conclude that the NRxx 
equations have the order of accuracy \(4M — 5)/3]. 

For boundary value problems, such discussion is only valid in the bulk. In the Knudsen 
layer, which is known to be of width O(Kn), we need to use X = x/Kn as the spatial 
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Figure 6: Shear and normal stress plots for the planar Couette flow with Kn = 0.5 



variable while investigating the accuracy of moment equations. In this case, if we consider 
a steady state problem, the small parameter no longer appears in the governing equations 
This means the order of magnitude for f a does not increase as \a\ increases, as 
has been found in [23] . 
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Figure 7: Density and temperature plots for the planar Couette flow with Kn = 1.0 

6.2 The validity of NR:e£c method for large Knudsen number and in the 
Knudsen layer 

As we have discussed above, there are two cases when the order of accuracy is not so 
meaningful for describing the accuracy of the NRxx method: 

1. In the case of Kn ~ 0(1), there are no "small parameters" in our concept. 

2. In the Knudsen layer, the orders of magnitude of moments do not increase as they 
behave in the bulk. 
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Figure 8: Shear and normal stress plots for the planar Couette flow with Kn = 1.0 



Nevertheless, we can still consider the NRxx method as a solver for Boltzmann equation 
with spectral expansion in the velocity space, and the method should be valid when f a 
decays sufficiently fast as |a| increases. Now we follow [23] and give the average absolute 
values of the moments with the same order for different Kn and different M in Figure [TTJ 
The result is based on the Couette flow problem in Section \5.2\ and the NRxx solution at 
x = —0.5 is used in the these plots. 

In these figures, we find that even when the Knudsen number is as large as 10, the 
magnitudes of moments still decay very fast. Thus, the NRxx method can still be consid- 
ered to be valid and efficient. Although the methodology of regularization which is based 
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;ure 9: Density, velocity and temperature plots for the planar Poiseuille flow 
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Figure 10: Stress and heat flux plots for the planar Poiseuille flow 
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on a small r is not valid any more, the regularization term (|2.1ip is simply a prediction 
of higher order moments. Such prediction differs from Grad equations' guess f a = 0, but 
also has a uniform expression for all Knudsen numbers. When M goes to infinity, the 
regularization term is expected to vanish since f a decays. On the other hand, this term 
smooths the profiles of the macroscopic variables, thus avoids the appearance of some 
unphysical phenomena such as subshocks (see [5]). This indicates the meaningfulness of 
regularization for practical use. 

7 Concluding remarks 

A uniform numerical scheme for coupling the NRxx method and the wall boundary 
conditions are developed in this paper, and the NRxx method is extended to apply the 
force term and predict correct Prandtl number by using the Shakhov collision model. To 
validate the proposed method, both steady and unsteady problems are simulated. We are 
currently working on applying the NRxx method to 2D problems. 
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Appendix 

A Some properties of Hermite polynomials 

The Hermite polynomials defined in (|2.8|) are a set of orthogonal polynomials over the 
domain (— oo, +oo). Their properties can be found in many mathematical handbooks such 
as pp. Some useful ones are listed below: 

1. Orthogonality: / He m ( ,)Jfc. W exp(-*V2) dx = 

2. Recursion relation: He n+ i(x) = xHe n (x) — nHe n -i(x); 

3. Differential relation: He' n (x) = nHe n _i{x). 

And the following equality can be derived from the last two relations: 

[He n (x) exp(-x 2 /2)]' = -He n+1 (x) exp(-x 2 /2). (A.l) 
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B Calculation of half-space integration 

The detailed calculation of I a ,p{8) (|4.9p will be presented. Using the definition of 
T~te,a(,v) (|2.7p . eq. (|4,9p can be rewritten as 



i*a ) = n 

where 



k=l 



(27T)- 1 / 2 ^-fffc / ,+ °° 

j — 2 / He ak (v k )Hep k (v k )ex.p 

a k [ Ji k 



\ v k\ 



h 



— oo, fc = 1, 3, 
0, jfe = 2. 



Applying the orthogonality of Hermite polynomials to (jB.ip , we have 

|2' 



(27T)^ 1 /2^^ 



,1 



He a2 (v 2 )Hep 2 (v 2 ) exp 



P2 



dt>2 



•<5, 



"2' JO 

Now it is obvious that (|4.1U|) holds if 

2 r+oo 

S(m, n) = — / He m (x)He n (x) exp(— x 2 /2) dx. 

\/2irm\ Jo 

Some simple cases can be directly worked out as 

i r+oo 

S(0,0) = -y=j exp(-x 2 /2)dx = 1/2, 

S(Q,n) = -= #e n (x)exp(-x 2 /2)dx = —=He n -i{$), n / 0, 
V2vr Jo V27T 

r+oo i 

S(m,0) = .— / ffe m (x)exp(-x 2 /2)dx = - ffe m -i(0), m^O. 

V27rm! Jo V27rm! 



(B.l) 



(B.2) 



(B.3) 



(B.4) 



(B.5) 



This agrees with the first three cases of (|4.1ip . For m/0 and n / 0, we use the differential 
relation of Hermite polynomials and get 



S^m, n) 



1 



2-7rm! Ae[o,+oo) 



He n (x) d[He m -i(x) exp(-x 2 /2)] 



1 



2irm\ 
1 



27rm! 



r+oo 

#e m _i(0)-ffe n (0) + n / iJe m _i(x)#e n _i(x) exp(-x 2 /2) dx 
Jo 

#e m _i (0)£Tc„(0) + n/m- S(m- l,n- 1). 

(B.6) 



This is the last case in (|4.1ip . 



C Expansion of the half-Maxwellian 



This section is devoted to the detailed calculation of p a defined in f|4.23j) . Due to 
(|4,26p . only (j4.24p and (|4.25p need to be evaluated. We first consider J s (x) with s ^ 1. 
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By applying the recursion relation of Hermite polynomials, we get 

J.{x) = j l 6 s * k -jL= exp (- 1 ^ ~ Xl 2 j ^.^(yJ-^-lJFe^CyJldi/ 

= — ^5-2(2:) + -J a -i(a;) 
s s 



2 

j! " 



/■ + °° 1 \ \V0y-x\ 2 \ ( 9 x r-\ TT 

(C.l) 

For the underlined term, we use integration by parts and the differential relation of Hermite 
polynomials, and get 

9 x 
J s (x) = --J s _2(x) + -J s -i(x) 

= -[(e w -e)j s ^ 2 (x) + xj s ^(x)]. 

s 

When s = or s = - 1, the integral ()4.27p can be directly worked out as (|4.30p since 
Heo(y) = 1 and He-\{y) = 0. 

The calculation of (|4.25p is almost the same as (|4.24|) . The only difference is that a 
boundary term will appear when integrating by parts. So the result becomes 



1 



1 e w 



X 



2 



J s (x) = - s [(9 W - 9)J s - 2 {x) + xJ s _ 1 {x)\ - ^y^0^Fe s _i(O)exp [~^w J, O 1 

(C.3) 

with initial conditions (|4.31|) . Define 



1 9 W s-i / x 2 \ 

H s (x) = -^—9—He s ^(0)exp \-^\ • (C.4) 

Then (j4.28j) and (|4.32p are natural. For s ^ 1 , the recursion relation of H s can be deduced 
as 



9 1 9 W s-3 ( x 2 \ 

(C.5) 
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Figure 11: Average values of {|/a|}| Q |=fc for k = 1 to 9. The results are based the NRxx 
solutions of the Couette flow. 
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